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Abstract: Time dependent perturbations of states in the holographic dual of a 3+1 
dimensional confining theory are considered. The perturbations are induced by vary¬ 
ing the coupling to the theory’s most relevant operator. The dual gravitational theory 
belongs to a class of Einstein-dilaton theories which exhibit a mass gap at zero tem¬ 
perature and a first order deconfining phase transition at finite temperature. The 
perturbation is realized in various thermal bulk solutions by specifying time dependent 
boundary conditions on the scalar, and we solve the fully backreacted Einstein-dilaton 
equations of motion subject to these boundary conditions. We compute the character¬ 
istic time scale of many thermalization processes, noting that in every case we examine, 
this time scale is determined by the imaginary part of the lowest lying quasi-normal 
mode of the final state black brane. We quantify the dependence of this final state on 
parameters of the quench, and construct a dynamical phase diagram. Eurther support 
for a universal scaling regime in the abrupt quench limit is provided. 
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1 Overview 

1.1 Holographic Thermalization 

To date, holographic theories of strongly interacting matter have been deformed, probed, 
perturbed and otherwise studied in an enormous number of theoretical experiments. 
The results of these investigations have provided many important lessons about the 
nature of strongly coupled field theories, even in applications that fall outside the well 
established examples of AdS/GET duality. While many earlier studies of holographic 
matter focused on the system’s linear response to the insertion of various operators, 
quite a bit of attention has recently turned towards understanding the full non-linear 
dynamics of these holographic field theories. 

Broadly, the aim of these applications is to uncover results about the processes 
by which strongly coupled gauge theories respond to arbitrary time dependent per¬ 
turbations. Often, as in the case of the ground-breaking early examples of numerical 
holography [1-4] , the gauge theory responds by “thermalizing” the perturbation. This 
means that after some characteristic time scale, typically set by the Hawking tempera¬ 
ture of a black hole in the dual gravity solution, the field theory arrives at a final static 
state in thermodynamic equilibrium. In this equilibrium state, all correlation functions 
assume their thermal values, which is to say the trace is taken with respect to the 
thermal probability distribution appropriate to the ensemble. In other examples [5-7], 
perturbations of the gauge theory have been found which apparently never thermalize. 
These “islands of stability” correspond to initial gravitational data that never leads to 
gravitational collapse and horizon formation. 

A conceptually simple means of perturbing a gauge theory is to turn on a time 
dependent source for some operator and evaluate the system’s response. Beyond the 
linear regime, such a source can be used to continuously drive the system, or to “quench” 
it. By the latter, one typically means that a coupling in the gauge theory is varied over 
a compact timescale f, where f is sometimes chosen to be small compared to other 
scales in the theory. 
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When quenched in this way, it is clear that the system’s response will depend on 
the properties of the operator to which the source couples as well as the parameters 
which specify the quench. In the weak field approach pioneered in [8], perturbative 
results have been obtained for the quench of a marginal operator in the case where the 
perturbation’s amplitude and characteristic timescale are both sufficiently small. This 
analysis allowed the authors to study the formation of black holes and black branes in 
asymptotically Anti-de Sitter spacetimes. 

Specifically, they constructed the limiting behavior of the dynamical phase diagram 
for the outcome of massless scalar collapse in global and Poincare patch AdS. Notably, 
they found that arbitrarily small perturbations of the Poincare patch solution always 
resulted in the formation of a black brane, whereas finite volume effects in the global 
solution prevented horizon formation under certain conditions. Interestingly, they also 
observe a crossover between final states with small and large black holes, as well as 
Choptuik behavior in the vicinity of the transition between final states with a horizon 
and those without. This scaling behavior strongly suggests the presence of a second 
order fine separating these phases. More recently, the weak field approach has been 
extended to provide interesting insights into thermalization in non-isotropic quenches 
[9], in finite density states of a gauge theory [10], and in a simple model of strongly 
coupled matter with a mass gap [11]. 

The holographic toolkit makes it conceptually straightforward to move beyond the 
weak field perturbative scheme as well as to generalize the quench to perturbations 
by relevant operators. Along the first direction, calculating beyond the reach of per¬ 
turbative methods clearly implies grappling with the full non-linearity of the Einstein 
equations. In general, this is an exercise in “numerical holography”, loosely defined as 
the set of all holographic calculations that require numerical techniques more involved 
than a call to Mathematica’s NDSOLVE. Instead, one typically appeals to any one of an 
assortment of numerical methods specifically tailored to the problem of gravitational 
in-fall. 

A recent example of holographic thermalization which connects perturbative weak 
field calculations to the full numerical evolution of a marginal perturbation appears in 
[7, 11]. In this setup, the authors consider the response of a confining gauge theory 
to time dependent perturbations by studying the dynamical evolution of a massless 
scalar quench in the AdS hardwall geometry. This gravitational theory is characterized 
by a length scale Zq = 1/K ok, which the geometry is artificially terminated. In turn, 
this length scale gives rise to a mass gap ~ A in the dual gauge theory, as well as a 
temperature above which a large black brane can appear in the bulk. Accordingly, the 
AdS hardwall serves as a crude model for a confining gauge theory with a deconfined 
phase dual to the large black brane solution. 
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Weak field calculations in the hardwall background indicated that the dynamical 
phase diagram ought to include a transition between perturbations which result in black 
brane formation and those that remain in a horizon-less scattering state. This should 
be contrasted with the situation in [8] in which such a distinction between final states 
was only possible when the dual field theory is defined on a sphere. Unfortunately, 
the weak field approach was unable to answer questions about the late time fate of the 
scattering solutions, which is of particular importance for the field theory interpretation. 
To resolve the late time behavior of the scattering solutions, a numerical method was 
introduced to solve the Einstein equations and evolve the system arbitrarily far forward 
in time. 

A particularly noteworthy result of the numerical investigation was the observa¬ 
tion that there exist perturbations such that the subsequent scattering solutions never 
undergo gravitational collapse. Since no horizon is formed in such a process, the im¬ 
plication is that the holographically dual strongly coupled matter does not achieve 
thermodynamic equilibrium at late times. This finding differs in important ways from 
the islands of stability identified in scalar collapse in global AdS. In particular, it demon¬ 
strates that a perturbation initiated by explicitly sourcing a marginal operator need 
not necessarily thermalize in the infinite volume gauge theory. To rephrase this result 
in a way that will be more directly applicable to our present work, the introduction of 
a mass gap in the dual gauge theory can strongly effect the thermalization time for a 
certain class of perturbations. 

A related, but logically distinct line of research involves the deformation of strongly 
coupled matter by time dependent perturbations of a relevant scalar operator. Numer¬ 
ical investigations into this situation appear in both systems with a finite density of 
charge carriers (for example [12]) and without [13-15]. We will focus on the latter 
scenario, in which an uncharged black brane solution is perturbed by varying the non- 
normalizable boundary mode of a massive bulk scalar in time. Holographically, the 
dual picture is that an initial state of the gauge theory in thermodynamic equilibrium 
is deformed by turning on a relevant operator over some timescale f. The strongly 
coupled matter generically responds by passing through a non-linear regime before set¬ 
tling once more into thermodynamic equilibrium, albeit in a different thermodynamic 
macrostate. 

One of the most interesting results to emerge from the numerics of [13-15] has 
been the appearance of a “universal fast quench regime” in which the change in energy 
density S after the quench scaled as a power law in the quench width f, 

^’fINAL — ^INIT ~ 2 A - d ' 
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This scaling is argued to manifest when the quench width is small compared to any 
other scales in the theory, such as the characteristic amplitude of the perturbation S. 
Moreover, the scaling was shown to be robust for many choices of the relevant operator’s 
conformal dimension A, and its appearance was further elucidated analytically in a 
recent series of papers [16-18]. 

Given the non-trivial dynamical phase structure and strong dependence of the 
thermalization time on perturbations in the hardwall model, as well as the appearance 
of universal regimes in the space of quench parameters, it is natural to wonder how these 
properties might manifest in a holographic dual to a gauge theory that is more similar 
to SU{3) Yang-Mills. In a sense we make more precise in the following subsection, 
this means we are interested in finding a natural way to “soften” the hardwall while 
retaining the mass gap ~ A that was responsible for the interesting dynamical features 
of the quench. The construction and subsequent perturbation of such a holographic 
model will be the focus of this work. 

1.2 Confining Gauge Theories and Holography 

Holographic models of confining gauge theories generically involve a gravitational the¬ 
ory whose solutions include those which explicitly break scale invariance in the bulk. 
This condition is necessary but not sufficient to allow for non-trivial temperature de¬ 
pendence of states in the dual gauge theory, as well as a mass gap at zero temperature. 
These models also include solutions whose metric breaks bulk conformal invariance, the 
most relevant of which will be branches of black hole solutions which are holographically 
dual to thermal states with non-zero stress-energy. The deconfinement transition will 
be a thermodynamic phase transition from solutions without a black hole, to solutions 
with one. 

In what follows, we will focus on a particular class of holographic models for con¬ 
fining gauge theories which are in the Einstein-Dilaton family. As this title implies, 
our model will involve a non-trivial scalar whose profile explicitly breaks conformal 
invariance in the radial direction, and can be tuned to produce gravitational solutions 
dual to states on either side of a deconfinement transition. Within this rather sim¬ 
ple class of models there exists a remarkable wealth of phenomenological possibilities, 
which include sharp deconfining transitions of first or second order, crossover behavior 
between phases, and theories driven from their UV fixed point by relevant operators of 
arbitrary dimension. 

When the gravitational theory is dimensionally reduced to five dimensions, all such 
bulk solutions as described above (with the exception of flows to a non-trivial IR CFT 
which are not relevant for holographic models of Yang-Mills) will have a naked singu¬ 
larity in the IR [19, 20]. The presence of this IR singularity is a signal that something 
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is missing in the effective holographic description. In known examples, typically the 
missing ingredients are an assortment of other bulk fields that will necessarily obtain 
vevs in the ground state solution. If such fields are reintroduced, the singularity is 
expected to be resolved. In [21], Gubser provides a simple criterion for a resolvable 
singularity which states that for a singularity to be resolvable it should be possible to 
construct a solution with an infinitesimal regular horizon surrounding the singularity. 
Such singularities are called “good” or resolvable singularities. 

Even in such “good” solutions, holographic calculations of correlators may be ill- 
defined in the presence of the naked IR singularity. This happens when in the neighbor¬ 
hood of the singularity the Sturm-Liouville problem that governs the bulk fluctuation 
has two normalizable solutions. In such a case an extra boundary condition is needed 
to determine the solution, and this is the signal that the correct solution is sensitive to 
the details of the singularity’s resolution. Such cases are called “holographically unre¬ 
liable” because without knowing the details of the singularity’s resolution holographic 
computations are inherently ambiguous. There are, however, numerous cases where 
the Sturm-Liouville problem near the naked singularity has only a single normalizable 
solution. In that case the solution is unique and does not depend on the resolution of 
the singularity. When this occurs, the singularity is said to be repulsive, [19, 20] and 
boundary theory correlation functions can be reliably calculated from holography in 
such backgrounds. 

In confining states of a gauge theory, the expectation value of Wilson loops exhibit 
area law scaling. Holographically, the Wilson loop has a natural definition as a minimal 
surface in the bulk affixed to the edges of a Wilson loop on the boundary. In [22], it 
was shown that holographic Wilson loops will give area law behavior providing that 
the string frame metric scale factor has a minimum, and that the scale factor at this 
minimum is non-zero. In Einstein-Dilaton theories, this requirement can be reformu¬ 
lated as a constraint on the IR behavior of the dilaton potential, as we will discuss in 
more detail below. 

In our model, the zero temperature, zero entropy^ ground state of the gauge theory 
will be dual to a bulk solution with a singularity in the far IR, but no event horizon. This 
singularity allows the bulk scalar to diverge in the IR, as well as the spatial metric scale 
factor to vanish.^ The latter is the feature responsible for vanishing entropy density, and 
is a familiar aspect of many candidate holographic ground states. As discussed above, 
the IR singularity can be thought of as a rather innocuous consequence of working in 

^Strictly speaking the entropy is 0(1) and is therefore due to one loop effects in the bulk. 

^It has been observed in [19, 20] that in Einstein-Dilaton theories which have genuine breaking of 
scale invariance, a discrete spectrum and a mass gap, the string frame metric becomes flat in the IR. 
Therefore, in the string frame the bulk “singularity” is due to a diverging dilaton. 
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a limit of non-critical string theory, providing the singularity satisfies several criteria. 
Our singularity is of this “good” type. 

As the temperature increases in the gauge theory, no new solutions are immediately 
available in the dual gravitational theory. Instead, the dual bulk solutions are given 
by the zero-temperature (confining) solution with its time coordinate compactified into 
the thermal circle. The temperature T associated to this solution is governed by the 
size of the thermal circle in the familiar way, 

t ^ ir where r = r +/9 with 15 = —, (1-2) 

and such a state is often called the “thermal gas” saddle point. As it lacks a black 
hole horizon, it is obvious that the thermal gas is characterized by thermodynamics 
subleading in the large (where Ac S> 1 is roughly the rank of the dual gauge theory) 
expansion inherent to our holographic setup. 

Increasing the temperature further, one eventually arrives at a special temperature 
T = Tq which marks the appearance of a new branch of solutions in the bulk theory. 
These solutions are of the black brane type, with a non-compact horizon of planar 
topology. As seen in the left plot of figure 3 for any T > Tq, there are two black hole 
solutions. One is in the left branch of the curve while the other is in the right branch 
of the curve. The black holes on the left branch are called large black holes as their 
horizon size is larger than any of those in the right-branch. They have positive specific 
heat and are therefore locally thermodynamically stable. On the contrary, the small 
black hole branch on the right has negative specific heat and the associated black holes 
are thermodynamically unstable. Qualitatively the diagram is similar to that of global 
AdS space, although here the space is flat and the volume infinite. In the case of global 
AdS, small black holes asymptote to Schwarzschild black holes in flat space in the limit 
of small horizon size (with T oo). Here also, small black holes in the limit of small 
horizon size have T —>■ oo but they are very different from Schwarzschild black holes, 
[20, 23, 24]. Finally, the unique black hole that is at T = To is very special, as here the 
specific heat diverges. 

The thermodynamics of such black holes is governed by the standard black-brane 
relations, which implies they have an entropy density proportional to L^/k^ where L 
is the characteristic length scale of the metric (for example the AdS radius) and k is 
the five-dimensional gravitational constant. A standard application of the holographic 
dictionary for known five-dimensional gravitational duals gives L^/k5 oc N^. The 
proportionality constant depends on the details of the duality under consideration. 
In holographic applications that do not descend from a known string theory solution, 
the proportionality constant is an undetermined parameter of the theory, but the 
scaling of various thermodynamic quantities is expected to be robust. 
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Figure 1. The phase diagram of our model. The thermal gas solution (grey) is the only 
solution below T = Tq. Above Tq the black hole branches appear (black plane), but the 
thermal gas is thermodynamically preferred. At T = Tg there is a first order deconfinement 
transition from the thermal gas to large black brane solutions (large black plane). Cartoon 
from [26]. 

As mentioned above, by studying the susceptibilities of the small black hole branch, 
for example the specihc heat, it is straightforward to show that these black brane so¬ 
lutions are thermodynamically unstable. This is a local statement, independent of the 
global thermodynamic (in)stability of these solutions. In light of the Correlated Stabil¬ 
ity Conjecture (CSC) [25], one might worry that these solutions are also dynamically 
unstable, which is to say their fluctuation spectrum may contain a mode in the upper- 
half complex frequency plane. Such an instability obviously leads to an exponential 
growth in time, and the conjecture posits that the existence of such a mode indicates 
a Gregory-Laflamme type instability towards a “lumpy” horizon. This dynamical in¬ 
stability will play no role in the present work as we will study s-wave perturbations of 
spatially homogeneous solutions. 

Above To, the thermal gas still provides the thermodynamically-dominant solution, 
as it has lower free energy than the black brane solutions at the same temperature. 
This remains true until T = Tg, at which point the model predicts a first order phase 
transition from the thermal gas to the black-brane solutions. Above Tg, the system is 
described by the large black brane solutions that are thermodynamically stable. This 
first order phase transition is a holographic realization of a deconfining transition to 
states with an entropy density that scales like and perimeter law behavior of the 
Wilson loop. The full phase structure of our model is illustrated in figure 1. 
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1.3 An Einstein-Dilaton Model of Thermalization in a Confining Theory 

Our goal in this work is to continue the investigation of the effects of confinement on 
thermalization processes in holographic gauge theories. We accomplish this by perturb¬ 
ing a state of a confining gauge theory with a relevant operator that we loosely associate 
with^ TrF^, and then studying the system’s response. From the bulk perspective, we 
are choosing a solution of our model and deforming it briefiy through time-dependent 
boundary conditions on the dilaton. The deformation will generically throw the sys¬ 
tem out of equilibrium, which will then undergo non-linear evolution governed by the 
Einstein-Dilaton equations. 

An important question is whether or not an arbitrary perturbation of a confined 
state in the strongly interacting gauge theory will thermalize. In this context, “to 
thermalize” means that the late time behavior of the state is in thermodynamic equi¬ 
librium, and all correlation functions are time independent and assume their thermal 
values. The gravitational formulation of this question is whether or not an arbitrary 
perturbation of the horizon-less thermal gas solutions of our model necessarily results in 
black brane formation. Put another way, is the exponentially diverging dilaton poten¬ 
tial of our model sufficiently similar to the hard wall of [7, 11] to encourage scattering 
solutions and prohibit horizon formation under certain conditions? 

At present, the answers to these questions remain beyond our grasp. In the absence 
of a horizon in the initial state, the diverging scalar in the IR leads to a number of 
challenging numerical obstacles. To address these issues, one should regulate the bulk 
solutions in the IR through the introduction of a numerical cut-off, and study the effects 
of the position of this cut-off on the computation’s results. A related approach, which 
we will adopt, forfeits the ability to comment on black hole formation while hopefully 
retaining qualitative effects of the confining dilaton potential. This scheme replaces the 
hard numerical cut-off with a small black brane horizon, which then acts as a natural 
IR regulator. 

In our model, this necessarily implies that our initial bulk solution will not be holo¬ 
graphically dual to the ground state of the gauge theory. At first glance it actually looks 
much worse, since these small black brane solutions aren’t the thermodynamically pre¬ 
ferred state at sufficiently low temperature, nor are they necessarily dynamically stable. 
Nonetheless, we argue that they provide a sensible starting point for our calculation, as 
in the limit of vanishing horizon size these black holes coincide with the ground state 

^In realistic bottom up holographic models of QCD, the scalar is dual to a marginally relevant 
operator [19, 20]. However, we do not expect qualitatively important changes in the IR if the UV 
dimension of that operator is three (as in the case studied here). The reason is that in both cases, in 
Einstein dilaton gravity the scalar becomes strongly relevant in the IR. 



of the system. 

The small black-brane solutions asymptotically approach the zero-temperature so¬ 
lution in the limit where the horizon recedes infinitely far from the bonndary into the 
IR. Thus, provided we accept the concession that a horizon has already been formed in 
the bulk, for sulRciently large separation between the UV boundary and the position of 
the horizon the gravitational solution shares many important features with the ground 
state solution. The most important of these are a scalar profile which is rapidly growing 
towards the IR, and a rapidly vanishing warp factor in the Einstein frame metric. We 
expect that these features will be responsible for most of the interesting dynamics in 
our system, such as the timescale with which thermalization takes place, and accord¬ 
ingly that the small black holes can provide us with a qualitative understanding of the 
conseqnences of these featnres. 

Moreover, the thermodynamic instability of these black branes and the possible 
presence of a dynamical instability in their fluctuation spectrum posited by the CSC 
is not expected to have much relevance to our discussion. Gregory-Laflamme type 
instabilities are related to a “clumping” of mass or charge at the horizon, and should 
thus involve finite wave-number perturbations in the directions transverse to the brane. 
We stndy translationally invariant pertnrbations in the non-critical (five dimensional) 
theory, and hence do not anticipate the possibility of exciting such instabilities. Indeed, 
we have not yet seen any direct evidence for a dynamical instability in the small black 
hole branch of our model. 

The main results of our study are a portion of the dynamical phase diagram, and 
the scaling properties of varions equilibration processes in onr model. The former is 
realized as a map between perturbation parameters and the final state achieved after 
equilibration. Our perturbation will be controlled by two parameters corresponding to 
the amplitnde of the pertnrbation and its duration. Thns, the dynamical phase diagram 
is a two dimensional plot with a curve marking the bonndary between perturbations 
that result in a small black hole and those that thermalize into large black holes. 

By focusing on the details of the final state, we are also able to make claims 
about the dependence of the characteristic thermalization time in the boundary gauge 
theory on the quench parameters. Not surprisingly, this timescale turns ont to be 
qnantified by the qnasi-normal mode spectrum of the final state solution. Moreover, 
by quantifying the relationship between the features of the quench and the change in 
energy density induced by the perturbation, we demonstrate explicitly that the simple 
scaling anticipated in (1.1) appears in our model as well. 

The remainder of this work is devoted to detailing the setup, implementation, and 
interpretation of time dependent quenches of a holographic confining gauge theory. In 
section 2 we introduce the specifics of the bulk theory we are interested in perturbing, 
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as well as the equations of motion that govern its behavior. The behavior of these 
solutions near the conformal boundary are of particular importance for the application 
of holographic methods, and thus we provide these near boundary solutions in some 
detail. 

The numerical construction and thermodynamic properties of the static, equilib¬ 
rium solutions of our model are discussed in section 3. This section primarily serves as 
an orientation to our holographic gauge theory. In it one hnds the dependence of the 
free energy and entropy density of our strongly coupled matter on temperature, as well 
as a comparison of the temperature dependence of the speed of sound in our model to 
that of SU (3) Yang-Mills. It is complemented by the results of section 4 which further 
characterize the static states in terms of their thermodynamic one-point functions. 

Section 5 details the computational approach we adopt to evolve our gravitational 
system in time. This approach is based on an assortment of well known numerical 
techniques that we carefully tune to accommodate the specifics of our model. Most 
notably, we explain how we handle the copious logarithmic fall-offs in the near boundary 
behavior of bulk fields. These fall offs are generic to gravitational theories in odd bulk 
dimensions, and present obstacles to the accurate determination of the normalizable 
and non-normalizable UV coefficients of the various bulk fields. This section will be 
primarily interesting to those who would like to numerically study dynamical quenches 
in related models. 

The output from our numerical method and the interpretation of this output is 
contained in sections 6 and 7. These sections constitute the primary results of our 
study. They include examples of the response our model to various classes of quench by 
a relevant scalar operator, as well as a dynamical phase diagram for the outcome of these 
quenches when performed in a particular initial state. We examine the dependence of 
the final equilibrium state on the quench parameters, and comment on the appearance 
of an anticipated universal scaling regime in the fast quench limit. The connection 
and applicability of these results to similar processes in other theories is discussed 
in section 8. Specifically, we comment on the implications of our calculations for the 
thermalization of probes in the strongly coupled matter produced in heavy ion collisions, 
as well as future directions one might wish to pursue. 

2 The Model 

2.1 The Action 

We will be interested in Einstein-Dilaton theories that are tuned to qualitatively repro¬ 
duce some important properties of QCD. These theories can be described by an action 
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of the form 

-5 = ^ / + V((f)^ Jg ( 2 . 1 ) 

where V{ip) is the dilaton potential which will govern the dynamics of the system, 
and /C is the Gibbons-Hawking-York term necessary to help render the variational 
problem well defined on the boundary. The model we study can be classified as a 
“bottom-up” description of a confining gauge theory, in the sense that V{ip) is not 
derived from a known supergravity theory. Instead, this potential is manufactured to 
satisfy certain criteria in the IR and/or UV limits of the gravitational theory so as to 
induce desirable features in the dual boundary theory. For example, these asymptotics 
determine whether or not there is a mass gap, while other asymptotics determine the 
scaling dimension of the dual scalar operator. These criteria have been elucidated and 
categorized in a series of papers beginning with [19, 20]. Importantly, the asymptotic 
behaviors of the scalar potential that we choose for the current work are present in 
known solutions of gauged supergravity. An example of this which is closely related to 
the model investigated in this work is described in [27]. As many qualitative features 
of holographic models are dominated by these asymptotic properties, it is reasonable 
to expect that predictions from our model may apply to a broad class of holographic 
systems, including those with a more esteemed string theory lineage. 

While there are a variety of scalar potentials whose distinct IR asymptotics lead to 
confining theories, it was found in [19, 20] that a particular IR behavior also produces 
linear radial trajectories for glueballs, as well as the structure of the Yang-Mills phase 
diagram just above the first order deconfinement transition. This behavior requires 

Vm ~ ( 2 . 2 ) 

which, apart from the subleading is the non-critical string theory dilaton potential 
in five dimensions in the Einstein frame. 

At high energies, the UV properties of the gauge theory depend on the nature of 
the scalar operator which perturbs it. In [19, 20] a marginal operator was used, and 
the potential chosen such that the UV fixed point was at (p — oo. The potential is 
of the form 

OO 

Uuv = ^ 

n=0 

so that A —> 0 where A in the UV is identified with the ’t Hooft coupling constant. This 
fixed point is extremely shallow as all finite derivatives of the potential vanish at the 
fixed point. 
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This potential is tuned to match the thermodynamics of SU (3) Yang-Mills, and the 
resulting theory, dubbed Improved Holographic QCD (IHQCD), is capable of quantita¬ 
tively reproducing many features of QCD across all energy scales. Clearly it would be 
desirable to study dynamical processes in IHQCD, beyond the linearized level. However, 
the presence of a marginally relevant operator in the UV creates serious computational 
difficulties. One consequence of the marginal deformation inherent to IHQCD is that 
the scalar A vanishes logarithmically near the UV boundary. This mild falloff is very 
difficult to handle numerically, as it requires retaining a very high level of numerical 
precision to correctly identify the leading and subleading coefficients governing the 
near boundary behavior of the solution. To circumvent this issue, we will sacrifice the 
quantitative match to Yang-Mills theory provided by IHQCD in favor of computational 
convenience. 

One can do this without departing drastically from the coarse features of Yang- 
Mills by following the approach advocated in [28]. In this approach, the marginal scalar 
operator dual to A is traded for a relevant operator with dimension not far from 4, which 
one hopes to roughly identify with a boundary operator of the form TrF^ where F is the 
Yang-Mills field strength. The fact that this operator is no longer marginal is meant to 
capture the anomalous dimension that the operator acquires after running some ways 
towards the IR, reminiscent of what happens in Yang-Mills. 

We will therefore assume that there is a regular UV fixed point at (/? = 0, without 
loss of generality, and the potential near this fixed point takes the standard form 

Uuv ~ + • • • , 9 ? ^ 0. (2.4) 

The various coefficients are fixed by the symmetries of the gauge theory and were 
shown in [19, 20, 29] to be in one to one correspondence to the coefficients of the 
holographic /3-function. The AdS scale L of the UV AdS space is determined by 
X 2 ]/(o) _ The relationship between the mass of the scalar and the conformal 
dimension of the dual gauge theory operator A can be made precise. In the standard 
quantization A is defined to be the larger of the two roots of 

= A(A-4). (2.5) 

8 

An example of a potential for which these properties are realized first appeared in 

[28]: 

12 (1-|-cosh |(/p — 

^^ ( 2 . 6 ) 

and a parameter set which produces a bulk theory dual to a confining gauge theory 
deformed by a dimension A = 3 operator is (a, b) = (1/500,10009/1500). This theory 
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is decidedly not real world YM. Unlike the more refined IHQCD models of [30], it fails 
to produce a very good match to the thermodynamics of real world SU (3) Yang-Mills 
(see section 3.1). Nevertheless, this holographic theory does exhibit a mass gap, the 
thermodynamic properties determined by its potential are qualitatively reminiscent of 
SU (3) glue, and it renders the bulk system relatively amenable to numerical investiga¬ 
tion. Accordingly, in what follows we will primarily be concerned with solutions to the 
equations of motion derived from (2.1) with the potential written in (2.6). 


2.2 Solutions 


Among the solutions to this model’s equations of motion are the planar geometries 
(with and without a horizon) of the form 

2 

ds^ = —Adv^ - -dvdz + S^da:^ and (p = (p{v, z) (2.7) 

Z 

where v is an ingoing null-coordinate, .2 is the radial direction, and x describe the planar 
E^. The various metric functions appearing in this ansatz, as well as the scalar, are 
taken to be functions of both v and ; 2 . The solutions satisfy the equations of motion 


0 = - 4 + 2(d+¥’)' + 3 ^ + 3 d+E k 

0 =E" + -s'+ A (, 3 “, 

0 = AV' + I cl+E T _ + A' + A", 


3z‘ 


E" 2= 3 
0=d5_E + i(d+^)=E + i2=d+E2l', 


( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 

( 2 . 12 ) 


where the prime (') denotes differentiation with respect to z, and d_|_ is a modified 
derivative defined by 

d+ = d,-jAd,, (2.13) 

which is the derivative along the outgoing null vector. 
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2.2.1 Boundary Analysis 

Near the UV boundary, the various metric functions and a bulk scalar with = —‘ijlT' 
can be expanded like 

A{v, z) = ^ K(n)+ anm{v) log"* z] 2 ;"-^ (2.14) 

n=0 m=l 

E(n, 2 ;) = [Sn(^’)+ ^nm(v) log"* z] z'^~^, (2.15) 

n=0 m=l 

^) = X] X ( 2 - 16 ) 

n=0 m=l 

and we adopt coordinates in which the asymptotically AdS^ boundary has ao(n) = 
So{v) = 1 and «om('*^) = <^om{v) = 0. In general, logarithmic terms are expected in even 
boundary theory dimensions, d, and when the scalar operator has dimension A such 
that A — d/2 is an integer. The scalar’s source, fo{v), remains a boundary condition 
to be implemented. 

Inserting these expansions into the Einstein equations and expanding near the 
boundary at 2 ; = 0 provides the asymptotic behaviors of the various fields. In what 
follows, we leave the potential largely unspecified, requiring only that it gives rise to a 
dual conformal gauge theory which is deformed by a parity invariant relevant operator 
of dimension A = 3. This constrains a few terms in an expansion of the dilaton 
potential about the AdS fixed point like 


V 

cp=0 


12 

L2’ 



for positive integers, n. Accordingly, the fields behave like this: 


(2.17) 


A 

E 


— ~ 2 ^ — -/q + a^z'^ — 0:4 z'^ log z + ..., 

= - + C “ g/o ^ + ^-^0 (^^Cfo — 4/0 j z"^ + S4^Z^ — a4^Z^ log z + ..., 
=/o + (^fo ~ C/o^ z'^ A f2Z^ — (f)2 z^ log z + . . . , 


(2.18) 

(2.19) 

( 2 . 20 ) 
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where 


04 - - g/o + g/oii) + /o ’ P'2') 

*4 = - 5/0/2 - ^fo + Y^/o (16C/0 + 3/„) + /„‘ , (2.22) 

04=i/„/. + /^(^-Y^V'<‘)). (2.23) 

/2 = -i/„ + /„^(ivW-^), (2.24) 

in which = d‘^V\^=Q. The functions a 4 , /o and /2 all depend on time, and 
are not determined by the asymptotic series expansion. The function C, is completely 
unhxed, and reflects residual reparametrization invariance in ; 2 . The coefficient 04 is 
related to the background energy density, and /o and /2 can be related to the source for 
and response of the dual scalar operator, respectively. The AdS radius L has been (and 
will continue to be) set to one, and can be reinstated in any formula by dimensional 
analysis. In order for the full set of Einstein equations to be consistently solved, the 
time derivative of 04 is subject to a constraint, 

=g {j 2 fo — fof 2 ^ + “^/oC (^/oC “ ^ (^4/0/0 — /o/o^ 

+ yC (/o' - /o/o) + /o/o' . (2.25) 


This constraint is related to the Ward identity governing the divergence of the stress 
tensor, which we revisit in section 4. 

From these expressions it is simple to work out the asymptotic behavior of the 
modified derivatives d+E and d+iy?. They are 


d+E—-(^- + (C^ —-fo ——fofo z + Cj: z'^ + ..., (2.26) 

d+9? == - 2/0 + ~ dC/o + |/o + ~ 2 ^^ z'^ + ..., (2.27) 


where (Es is the exhausting constant 


& = Ya4 + g/ 0/2 - gCVo + 4/0 (llC/o 


-/o/o - /n f — +- I • (2.28) 

216 1152 J ^ ’ 
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For later convenience, we write these asymptotic expansions in the form 

d+E(n, z) = '^ [{Ds)niv)+ {Da)nm{v) log"" (2.29) 

n=0 m=l 

d+^{v,z) = E muv)+ {D(j))nm{v) log"" z] z". (2.30) 

n=0 m=l 

Solutions with non-zero /o correspond to adding to the boundary Lagrangian a 
term like 

SC = fo O, where [/o] = (2.31) 

for mass scale A. For the model of interest (2.6), we have A = 3 and thus the source 
for the relevant perturbation has dimensions of energy. Regularity of the solutions in 
the IR implies one constraint between the UV coefficients fo, f 2 and a 4 , and thus the 
UV deformations are specified by a single dimensionless quantity which we can take 
to be {Ttt)/ff. In other words, we anticipate a one parameter family of black hole 
solutions characterized by their energy density in units of fo. Often times it will be 
more instructive to parametrize various solutions by the value their scalar obtains in 
the IR, a method introduced in the following section. 

3 The Initial State 

We eventually wish to thermalize the scalar perturbation in an initially static state 
of the gauge theory dual to the bulk theory described by (2.1) and (2.6), at finite 
temperature. These states are constructed by restricting the background functions 
in the ansatz to vary only radially, that is removing n-dependent terms in (2.8-2.12), 
and then integrating the resulting system of ordinary differential equations to obtain 
solutions for A{z), E(z) and (p{z). These bulk solutions then provide the starting point 
for the subsequent deformation and evolution. 

To construct numerical solutions to the equations of motion (2.8-2.12) describing 
the initial state, we exploit the fact that we desire a regular horizon at z = zh and thus 
expect that the background functions can be expanded around the horizon like 

F{z) ^ Fo + Fi{zh - z) + F 2 {zh - zf + ..., (3.1) 

where F is any of {A, E,(p}. The coefficients Fi can be partialiy fixed by appeafing 
to symmetries of the underfying gravity theory. For example, Aq = 0 by assumption 
of a regular horizon, and Sq = 1 can be enforced by rescalings of the spatial coordi¬ 
nates. As mentioned previously, the metric ansatz (2.7) does not fully fix the gauge 
freedoms of the gravitational theory, and this is reflected in horizon data as the ability 
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to choose Ai arbitrarily. In practice, we have found it convenient to fix this residual 
reparametrization invariance by requiring that the boundary coefficient ( of the previ¬ 
ous section vanishes. The value of the scalar at the horizon, (/?o = thus parametrizes 
the available solutions in this model. All the higher order coefficients are fixed in terms 
of these first few coefficients, and accordingly one can use the expansion to sufficiently 
high order to provide IR data at a radial location near the horizon with which to seed 
a numerical routine. 


3.1 Thermodynamics 

The solutions to the equations of motion are asymptotically AdSs near the boundary 
at z = 0, with a scalar that vanishes linearly. However, as our integration strategy 
relies on fixing the values of various parameters at the horizon, it will often happen 
that the numerical solutions obtained are characterized by different energy scales, /o 
and appear in different coordinates. More explicitly, static solutions near the boundary 
generically behave as 

=-(? + ■■■) - |d«U- + (If 

if =foz + ... (3.3) 

for constants T,p and /q. As the aim is to compare states of different temperature 
in the same theory, we should arrange that the solutions take a canonical form at 
the boundary, with identical metric and energy scale. This can be accomplished by 
performing a coordinate transformation to a new radial coordinate in which /o = 1, as 
well as simple rescalings of the other coordinates. To wit, the transformation 



V =/o U, 

X =foT,F X, 

Z =foZ, 


ensures a near boundary solution of the form 


ds2 




(f=Z+..., 


(3.4) 

(3.5) 

(3.6) 


(3.7) 

(3.8) 


which is then suitable for comparing thermodynamic properties between solutions. 
Following the standard prescriptions for black brane thermodynamics, temperatures 
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Figure 2. The free energy density in units of the critical temperature Tc as a function of 
T/Tc for the model described by (2.6). 


and entropies can be readily extracted from the near horizon geometry. In the “canon¬ 
ical” coordinates of (3. 4-3. 6 ), the temperature is easily computed from the surface 
gravity k at the horizon since 


r = 


k 

2it 


where 


zh 


(3.9) 


and is a unit Killing vector, timelike outside the Killing horizon dX z — Zh> Meanwhile, 
the entropy density is proportional to the area of the horizon^, so that 


T 


—zlA 

47r ^ 


zh 


and 



(3.10) 


where s is the entropy density and 7 is the determinant of the spatial part of the metric 
at the horizon. For the solutions we consider, inserting the near horizon form of the 
canonical metric into (3.10) gives 


T = 


4 vr/o 






and 



(3.11) 


As different solutions are fully characterized by the value of their scalar at the 
horizon, it is often useful to think of their thermodynamic properties as functions of 
\h = exp(/ 7 j^. In this spirit, once the temperature and entropy of a given state can be 
reliably computed, one can measure the free energy of the solution from the integrated 
first law: 

H^h)= (3.12) 

J\h dAif 

^For the sake of simplicity, we will often refer to the “area” of a horizon when we more precisely 
mean the “area density”. All of the black brane solutions we discuss have planar horizons and thus, 
strictly speaking, infinite horizon area. 
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Figure 3. Plots of the temperature scaled by the critical temperature as a fuuctiou of Xh/^c 
(left) aud the eutropy deusity scaled by the third power of the temperature as a fuuctiou of 
T/Tc (right). The rightmost plot becomes “dotted” as oue passes through the phase trausitiou 
by loweriug the temperature from above. This is meaut to iudicate that iu the field theory, this 
low temperature phase is goverued by the thermal gas solutious, whose eutropy is subleadiug 
iu the uumber of colors Nc. 
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Figure 4. Plot of the temperature dependence of the speed of sound. The bine curve is the 
gravity prediction for the model described by ( 2 . 6 ), while the red sqnares are lattice data for 
SU{3) gauge theory from [31]. The blue dashed line indicates the conformal value = 1/3. 

For the model considered here, the free energy density is plotted as a function of T in 
figure 2. 

From these plots, one learns that there exists a critical \h = Ac, located where 
the free energy changes sign, at which there is a first order phase transition from the 
thermal gas to the black brane phase. One also notes the presence of two black hole 
branches in the plot of as a function of T. Since the specific heat is given by 

C, = -T-^, (3.13) 

wherein T is manifestly positive, one finds that the “upper” branch of black hole solu- 
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tions, with positive curvature, have < 0 and are thus thermodynamically unstable. 
Accordingly we will anoint this the “small black hole” branch in analogy to similar 
solutions in global AdS. In hgure 3 the temperature is plotted as a function of \h, and 
the entropy density is plotted against the temperature. From the former one notes the 
existence of a minimum temperature, which we denote Tq, while from the latter one 
finds that this model is characterized by a large jump in the entropy density at the first 
order phase transition. Moreover, as the squared speed of sound is just 


2 ^ dlogT 
® d log s ’ 


(3.14) 


this model attains the high-T conformal value of = 1/3 more rapidly than anticipated 
by SU (3) glue on the lattice. The results are shown in figure 4. 

In what follows, we will focus our attention on finite temperature initial states in 
thermal equilibrium. This is a subset of thermalization processes which excludes the 
possibility of commenting on a variety of interesting questions related to the conditions 
required for the formation of a black hole in this theory. Included in this subset, 
however, are processes which probe the non-trivial gravitational phase structure of the 
bulk theory. Generically, one expects that this phase structure will result in an interplay 
between physics on the large and small black hole branches, which in turn might better 
inform our understanding of thermalization in confining gauge theories. 


4 More Boundary Theory Observables 

The dual gauge theory data is encoded in correlation functions of the boundary theory 
operators. In this work, we will primarily be interested in the one-point functions 
of the stress-energy tensor and the dimension three operator O. Holographically, 
the values of these one point functions are related to the boundary coefficients of the 
normalizable modes of the metric and bulk scalar, a 4 and /2 respectively. 

4.1 Renormalized One-Point Functions 

The precise relationship, however, requires a careful analysis of the near boundary 
onshell bulk action, which generically has both power law and logarithmic divergences 
at z = 0. These divergences can be regularized and renormalized following the standard 
dogma of Holographic Renormalization [32, 33], the end result being a (possibly scheme 
dependent) identification of bulk falloffs with boundary theory correlation functions. 
This procedure is by now a familiar aspect of many holographic calculations, and 
accordingly the details will be left to appendix A. 
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The main results of this analysis are the set of local counterterms required to 
regulate the on shell action, and the one point functions of the stress energy tensor 
and the scalar operator as functions of the bulk boundary data. In the present system, 
with flat boundary metric and (p dual to a dimension three operator, the relevant 
counterterms turn out to be 


5', 


ct 


= 6 + + loge^F4<^^ - +A[j,^] 


(4.1) 


where 7 is the pull back of the metric to the radial cutoff at e, .A is a set of finite 
counterterms which define the renormalization scheme, and F 4 = 16/27 — j2^. 

Absent supersymmetry, or another guiding principle with which to fix the coefficients 
of the finite counterterms, we shall henceforth adopt a holographic minimal subtraction 
scheme and simply ignore them. The corresponding one point functions are then given 
by 



1 2/ 1-\ 4- 1 /14 f/(T\ 

^ - 2^4 + g [2f2 + -foj /o + ^/o + g/o , (4.4) 


where we have used (A.25) to express the correlation functions in terms of the co¬ 
efficients in (2.20). As expected, these correlation functions are constrained by the 
presence of Ward identities. In terms of these near boundary expansion coefficients, 
the conformal Ward identity reads 

(r,) = ? ( 4/2 - /„) /„ + = MO) + ?/„" - (4.5) 

while the divergence of the stress tensor yields 


V‘(T„} = MO), 


(4.6) 


which illustrates the non-conservation of the system’s energy in the presence of time 
dependent sources. 

To better interpret the values of these various one-point functions, it is convenient 
to define subtracted correlators which effectively measure deviations of the energy, 
pressure, and/or scalar expectation value from the static zero temperature solution 
obtained in the limit ipn 00 . One way to do this is to modify directly the renormal¬ 
ized on shell action by addition of appropriate finite local counterterms. For example. 
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Figure 5. “Hatted” one point functions characterizing the static states of the boundary 
gauge theory. The plots are given in units of /o with k? — 1. The dotted (left) and the 
dashed (right) lines correspond to T — Tc and T — Tq, respectively. 


a non-zero A = introduces a parameter which can be tuned such that the 

energy density of the zero temperature solution vanishes. Since the counterterms are 
properties of the theory and not a specihc solution, this subtraction will be manifest in 
all one-point functions computed in the boundary theory. 

A related, and perhaps more practical prescription is to simply define new correla¬ 
tors with the limiting ipn oo values explicitly removed. Accordingly, in the examples 
that follow we often construct “hatted” one-point functions defined such that 

(Cl) = (Q) — {^)ipjj=cx> (4.7) 

for any boundary theory operator fl. In figure 5 we plot the hatted correlators as 
functions of the value the scalar obtains in the IR for future reference. 

The numerical procedure we adopt is initialized by the boundary coefficients 04 , /2 
and /o, which are in turn extracted from numerically generated initial states as de¬ 
scribed in the previous section. For this reason, the evolution’s stability depends cru¬ 
cially on the accuracy of these values. It is thus reassuring that we find excellent 
agreement between the free energy density computed from ( 3 . 12 ), and from boundary 
data using T = —p = —{T^x)- 

We also compute the dependence of the energy density (T)t) on the system’s entropy 
density s. The former is calculated from the UV boundary coefficients extracted from 
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Figure 6. The energy density {Tu) as a function of entropy s, in units of /o and with /^ = 1. 
The asymptotic behaviors are {Tu) oc and {Tu) oc Ins in the limits of very large and 
very small black holes, and a fit for the former is plotted with a red dotted line. The green 
and magenta dots mark the locations of the first order phase transition at T = Tc and the 
division between small and large black holes at T = Tq, respectively. 

our numerical solutions, while the latter is computed from horizon data. The result 
appears in figure 6, in units of /o and k = 1 . At large energy density, our model correctly 
reproduces the expectation for a conformal theory, {Tu) oc 5^/^. As the energy density 
decreases, the dimensionful source’s explicit breaking of the UV theory’s conformal 
invariance becomes increasingly important. In the extremal limit, the small black hole 
branch is characterized by a logarithmic dependence of the energy density on entropy 
density, {Tu) oc s-\/— Ins. From the point of view of the boundary gauge theory, the 
latter behavior does not manifest as the thermodynamically preferred solution is the 
thermal gas whose entropy density is subleading in Nc. 

5 Numerical Strategy 

The primary goal of the present work is to discover how an initial state described by 
a solution from section 3.1 responds to time dependent perturbation by a dimension 
three operator. Holographically, we implement this perturbation by varying the UV 
boundary condition on the scalar in time. For example, a source whose profile is given 
by 

fo{v) = fo - Sfoe~^ (5.1) 

corresponds to perturbing the initial state (with coupling to O fixed by /o = /o(—oo)) 
by a gaussian deformation with amplitude —(5/o and variance centered at n = 0. 
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Accordingly, the boundary theory experiment we have in mind is dialing down the 
coupling to the relevant operator by an amount Sfo and for a time r. 

In the remainder of this section, a recipe for computing the time evolution of an 
initial state in the presence of time dependent scalar source fo{v) is explained from 
both generic and practical viewpoints. 

5.1 Integration Routine 

The equations of motion (2.8-2.12) have been arranged into a “nested” structure con¬ 
venient for numerical study. We employ and consequently review the method detailed 
in [34], which exploits this nested structure to reduce the partial differential equations 
to a sequence of ordinary differential equations that can be solved time step by time 
step. 

Before addressing the task of solving the Einstein equations describing this system, 
it is important to understand the boundary information needed to fully specify a given 
solution. For this, one may turn to the near boundary behavior of the homogeneous 
Einstein equations. For example, in the UV equation (2.10) is 

S" + -S' + ^/2S = 0, (5.2) 

y 

which has the general solution 

^ ^ 1 2foz I . 2foz 

L = Ci-cos—-|-C2-sm—-—. (5.3) 

Z 6 Z 6 

Thus, near the boundary linearly independent solutions behave like z~^ and Com¬ 
paring to (2.19) it is clear that E can be fully specified by the first two terms in the 
UV expansion—the AdS boundary condition and the radial reparametrization artifact 
C,. Similarly, from the unsourced (2.9) one finds near the boundary 

(d+S)'--d+S = 0, (5.4) 

Z 

which has the solution d+E ~ z^. From (2.26) it is clear that one must specify the value 
of Cs in order to determine d+E, which in turn depends on /o, /2, and 04. Finally, 
from the unsourced (2.8) in the UV one solves 

(d+(p)'- ^ d+(p = 0 (5.5) 

with d+(y? ~ z^!"^. From (2.27), evidently the desired solution is the one in which the 
coefficient of z^!"^ vanishes. In principle, the same sort of analysis can be applied to 
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( 2 . 11 ) to explore the boundary behavior of A. The homogeneous equation is readily 
solved by A ~ z~^ + z^, which (2.18) shows are fixed by ( and ( on the given time 
slice. Since ( parameterizes a residual gauge degree of freedom, one may adopt a gauge 
where C = 0- Alternatively, it may be desirable to determine ((v) dynamically. In this 
case one may use the value of A on the apparent horizon instead of ( as an integration 
constant. In practice this can be accomplished by immobilizing the location of the 
horizon (e.g. requiring Zff = 1 always). The condition for the location of an apparent 
horizon in these backgrounds (see Appendix C) is simply 


d+S = 0 , 

^ I ZH ’ 

and the horizon stationarity equation thus requires 


A(zh) 


16 (d+(p)^ 
3 V 


(5.6) 


(5.7) 


The upshot of the preceding analysis is that solving the Einstein equations requires 
knowledge of 04 , ( and (p(vo, z )—from which /2 can be extracted—on the initial time 
slice, together with a choice of the forcing function fo(v). 

After obtaining this initial data, one is well poised to evolve the system. The 
procedure is as follows: 


1 . From (p(vo,z), the AdS boundary condition So(^o) = and the value of ((vq), 
equation ( 2 . 10 ) can be integrated to obtain E(uo, 2 :). 

2. With this knowledge, and the value of Cs(uo) —which depends on 04(^0), /o(n) 
and f 2 {vo )—equation (2.9) can be integrated for d+E(uo, 2 :). 

3. Then equation (2.8) can be solved with the boundary condition that d +99 has no 
fall off like at the boundary, resulting in the knowledge of d+(p(uo,-s). 

4. Finally, we turn to the second order ODE given by equation (2.11). As indicated 
above, this equation can be solved given the values of C(^o) and either C{vq) or 
Ah{vo). 


Already this routine is enough to evolve the fields E and ip, as well as ( in time. 
More precisely, for any field F, knowledge of d+F permits one to write 

F = d+F H- AF' and thus F{vo + An) ~ F(uo) + F{vo)Av. (5.8) 
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In practice, it is not necessary to evolve S explicitly, as its profile on the next time slice 
will be constructed when step 1 repeats. The value of f 2 {vo + An) can be computed 
by extracting f 2 {vo) from the boundary behavior of (p(no, z) and integrating, or it can 
be extracted directly from the boundary behavior of (p(no + An). In practice we have 
found the latter numerically favorable. 

At this juncture, all that is needed to update the routine on the next time step is 
ah (no). There are options for computing this. One method is to solve equation ( 2 . 12 ), 
rearrange the d+ derivative to obtain 9^d+S at no, and then study the UV behavior to 
extract a 4 (no). A better way is to use (2.25) and directly evolve a 4 (no) to the next time 
slice. In this way one arrives at time n = no + An with updated values of / 2 , 04 , and 
(p{vo + Av,z), which is all the information required to begin the integration routine 
anew. 

5.2 A practical method 

The radial integration of the equations of motion can in principle be performed by any 
of a number of well worn techniques, including pseudo spectral and finite difference 
methods. In practice, we found that the logarithmic terms in the near boundary be¬ 
haviors of the fields rendered pseudo spectral methods unstable unless sufficient care 
was exercised to explicitly remove these terms (an observation also made by the authors 
of [15]). While one can typically ameliorate this issue through suitable field redefini¬ 
tions, we have found it more convenient to employ a finite difference discretization 
method and integrate from the boundary to the IR. This method is similar in spirit to 
the one that appears in [35]. 

The response of the boundary gauge theory to these time dependent perturbations 
is encoded in the time evolution of /2 and a 4 , which show up in higher orders of the 
boundary series solutions (2.18-2.20). Therefore, to determine them accurately, we 
work with redefined fields where /2 and 04 appear in the leading order at the boundary 
expansions by subtracting source contributions in lower orders. To surmount some of 
the challenges posed by the logarithmic behavior near the boundary, we also subtract 
the first few logarithmic terms. Moreover, we rescale by appropriate factors of .2 so that 
the redefined fields do not vanish at the boundary. We thus work with the subtracted 
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fields defined by 


3 6 


z'^ A = A — ^ aniz^ ^ log — 062 ^^ log^ ; 2 , 


(5.9) 



2^ E = E - S„2" * - ^ fT„i2” ^ log Z - 0622^ log^ 2, 


( 6 . 10 ) 


n=0 n=4 

1 4 


Z^tp^tp - y^/„2”+' - y^'(>„l2"'^' log2 - 1)1422® log^ 2, 


( 5 . 11 ) 


n=0 n=2 

3 6 


= d+E - Y,^Ds)nZ^-^ - Y.iDaUz'^-Hogz - {Da)^ 2 Z^\og^ z, (5.12) 



Dip = d+p - '^{Df)nz'^ - ^(D0)„iz” log z - (E>0)42Z^ log^ z. (5.13) 


In our choice, the redefined fields are (7^ over the radial domain. Substituting these 
expressions into the equations of motion, we obtain evolution equations of the redefined 
fields, which are solved with the strategy described in section 5.1. 

The equations of motion are discretized by finite difference grids in the z and v 
coordinates whose sizes are Az and An, and we arrange our numerical computation so 
that second order accuracy is supposed to be achieved. We use a fourth-order Runge- 
Kutta method in z-integration, and integrate from the boundary toward the interior. 
The z-derivatives of the fields which are not evolved in each differential equation are 
computed by a central finite difference scheme. We adopt a gauge where ((^(u) = 0 is 
fixed for all times. This generically implies that the location of the black hole horizon 
can vary throughout the evolution. We monitor the location of the apparent horizon on 
each time step and terminate the integration slightly inside it, which is also inside the 
event horizon. To update (p to the next time step, we use an upwind difference scheme 
for the advection term in (2.13). The grid sizes are chosen such that the Courant- 
Friedrichs-Lewy condition is satisfied, and we use An < Az. For u-evolution, we use a 
modified Euler’s method, where to update 04 we integrate the constraint (2.25) with a 
fourth order Adams-Bashforth formula,^ 


55a4(u,j) - 59a4(u„_i) + 24a4(u,j_2) - 9a4(u,j_3) 

24 


«4(W+l) = «4(W) + 


An. (5.14) 


Once the final static configuration is reached, when the apparent and event horizons 


coincide, we solve (C.8) backward in time to compute the evolution of the event horizon 


®For first few steps, (5.8) is used. 
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Figure 7. The area density of small black holes compared with that at the phase transition, 
Ahc- The red dot marks the location in our space of solutions of the smallest black hole we 
perturb in this study. The dashed line indicates the location of the small and large black 
hole transition at T = Tq. As these are static black brane solutions, the apparent and event 
horizons coincide. 

^eh(^)- We can then calculate the area density of the event horizon on constant u-slices 
as 


^eh('^) — S('u,zeh)^- (5.15) 

The monotonic increase of the event horizon is utilized as a consistency check of our 
numerical computation. The area density of the apparent horizon is computed analo- 
gously, as AAuiv) = T,{v,ZAHf- 

6 Thermalization Examples 

In this section we provide examples of the thermalization processes produced by our 
numerical computation. As the scalar source of the static states is nonzero, quenches 
are not symmetric between positive and negative h/o, and we consider quenches with 
Sfo > 0 for simplicity. Several tests suggest that broad stroke, qualitative features of 
the resulting time evolution seem not to be significantly altered by the choice of this 
sign.® Throughout this section and the next we will work in units of /o and henceforth 
set = 1. The stability and accuracy of initial data is discussed in Appendix B. 

An important question in the interpretation of these results is “how far” our initial 
states are along a given black hole branch. A practical metric for quantifying this 
distance is the area density of the black hole horizon. The zero temperature solution, 
which can be thought of as the Xh ^ oo limit of the small black hole branch, has 

®The primary differences occur for very large, very slow quenches but we will primarily be interested 
in more modest perturbations in what follows. 
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Figure 8. Typical behavior of the magnitude of the late time deviation of {Txx) from its final 
value. This particular ring-down corresponds to the late time behavior of the perturbation 
shown in figure 12. The red line is a fit to the exponential decay provided by 5{Txx) ~ 
0.65 

vanishing horizon area, and the horizon area grows monotonically as the scalar vanishes 
at ^ 1. This further justifies our “large” and “small” naming conventions for the 
black hole branches. In the studies which follow, the smallest black hole we perturb 
has Xh/K = 3.23. In figure 7, the area of the black hole horizons as a function of A 
is given, and this smallest black hole is indicated by a red dot. Evidently, the horizon 
area of Xh/Xc = 3.23 in this case is about 1/5 of that at the first order phase transition. 

Generically, a time dependent perturbation of a black hole will take a static initial 
state through a non-linear regime initiated by the details of the quench profile, followed 
by a linear regime governed by the “ring-down” to the final steady state configuration 
(there may also be late time power law tails in some situations, but we will not be 
concerned with these here). The ring-down is fully determined by the quasi-normal 
modes of the final state black hole, and is dominated by the mode closest to the real 
axis, oji- In turn, the quasi-normal modes characterize the linear response of the final 
state to small perturbations in any of several available channels. In the present case, 
where the perturbations preserve the homogeneity of the spatial M^, the gauge invariant 
perturbations organize themselves into representations of SO{3) [36, 37] transforming 
as the transverse-traceless (spin-2), vector, or scalar. 

In figure 8 the late time ring-down of one particular time dependent perturbation 
is shown on a logarithmic scale. The figure clearly indicates the presence of an excited 
mode of the form 


8{Txx) ~ ReZie with tui = tu* — iV (6.1) 

for some real constants Zi, and T > 0. Although figure 8 illustrates the general 
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Figure 9. The temperature dependence of the decay width T for the lowest lying scalar quasi¬ 
normal mode in several states of our theory. The blue circles are large black branes whose 
temperature is an integer multiple of Tg. The orange squares correspond to the minimum 
temperature black brane (top) and the smallest black hole we perturb in our study (bottom). 
The ratio T/ttT approaches 1.75953 (the dashed line) at high temperatures, which coincides 
with the expected value for perturbations of AdSs Schwarzschild by a dimension 3 scalar 
operator [38]. 


late time features of any perturbation in our study, it is important to note that the 
quasi-normal mode spectrum, and hence the particular value of tui, is a property of the 
final state achieved after the quench. In other words, one should keep in mind that ui 
varies as one traverses the static state solution space so that uji = C(Ji(Aj^). In figure 9 
we quantify this by plotting F as a function of temperature for several thermal states 
of our holographic theory. As expected, linear scaling of F with temperature appears 
in the conformal (small scalar) limit, and deviations from this behavior are already 
evident at the first order phase transition where T = Tc. The qualitative properties of 
this plot are anticipated by the thermodynamic features of our model shown in figure 
4 and the right plot of figure 3 which display a similar approach to conformal behavior 
above Tc. 

To construct a meaningful measure of the “thermalization time” in our quench 
process, it is useful to imagine partitioning the system’s response to a time-dependent 
perturbation into several distinct timescales. The gaussian perturbations that we de¬ 
fine in 5.1 are controlled by two dimensionless parameters, S = (5/o//o and f = r/o. 
Thus, the timescale f is roughly the amount of time it takes to drive the system out of 
equilibrium. Once the system is no longer in thermodynamic equilibrium, it may return 
to a (generically different) static thermal state by passing through any number of addi¬ 
tional dynamical regimes, each with its own characteristic timescale, 71- For example, 
one might plausibly imagine a regime with timescale 7s dominated by scattering of the 


30 




perturbation in the bulk from the rapidly increasing scalar potential, or a regime whose 
characteristic timescale 7rd is controlled entirely by the low-lying quasi-normal modes 
of the final state black brane. Of course the system could respond in more complicated 
ways, and the distinction between regimes may not necessarily be crisp. In any case, 
we define the thermalization time as the sum of these characteristic timescales: 

Ttherm = %■ (6.2) 

i 

An important attribute of this definition is the fact that f is excluded from the sum. 
This feature ensures that the thermalization time we measure is characteristic of the 
system’s response to a given perturbation, and does not explicitly depend on the time 
it takes us to prepare the non-equilibrium state. 

In practice, we find that our quenches are well described by a very rapid transition 
from the regime driven by the quench to the ring down characterized by linear response. 
Put another way, in (6.2) only one term, 7rd, appears in the sum. We will thus be 
particularly interested in the decay width P, as it provides a natural time scale for all 
of the equilibration processes we examine. In other words, throughout the rest of this 
work we will always find that 

Ttherm ~ ~ • (6.3) 

r((5,r) 

Given the preceding discussion, there are two natural questions we would like to address: 

1. For a given initial state, what is the dependence of the thermalization time on h 
and f ? 

2. Why does this thermalization time appear to be dominated by 7rd in our con¬ 
fining gauge theory? 

Our answers to these questions motivate the remainder of this work. The answer to the 
first question begins in the following section, while we postpone addressing the second 
until the discussion in section 8. 

6.1 Thermalization Between Branches 

The form of the Ward identity (4.6) encourages the intuitive expectation that a large, 
fast quench will result in the most significant increase in the system’s energy. In 
figure 10, examples of the evolution of the boundary operators in such a scenario are 
shown. These computations are performed with a width f = 0.168 and amplitudes 5 = 
0.5 and 1. These perturbations basically dominate the dynamics, and can generically be 
tuned to result in a final configuration characterized by a large black hole independent 
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Figure 10. Large amplitude quench. The blue and purple lines correspond to 5 = 0.5 and 
1, respectively. In the case of the larger amplitude quench, it is interesting to note that the 
energy density appears to be driven below the ground state energy density (i.e. negative) in 
the first moments of the quench. 

of the particulars of the initial state. Indeed, we have performed analogous quenches 
in many different initial states, all of which show qualitatively similar time evolution 
during and after the quench. In line with our previous discussion, the plots of {O) 
and (Tjj) show late time behavior that is well fit by an exponential decay related to 
the ring-down of the black hole. As expected, we find that the lowest lying quasi¬ 
normal mode for these large black hole final states retreats from the real axis as the 
final state energy density increases. Accordingly, we observe that these large amplitude 
perturbations imply rapid thermalization in the dual field theory. 

Once the entire evolution between static initial and final states is known, the time 
evolution of the apparent and event horizons can be computed. This provides a useful 
view of the equilibration process, as well as a consistency check of the time evolution. 
The area of the event horizon monotonically increases as shown in figure 11, and its 
location should always lie outside (closer to the boundary than) that of the apparent 
horizon. We verify explicitly that both of these statements hold in each process we 
study. Additionally, we note that in all perturbations we evolve, the area of the apparent 
horizon increases monotonically as well. 
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Figure 11. Example of the time evolution of the apparent (red) and event (blue) horizons. 
The event horizon coincides with the apparent horizon when the bulk solution is static, at 
V Too. 



V V 
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Figure 12. A “small” amplitude quench with = 0.02 and f = 0.168. 

6.2 Thermalization Along a Branch 

When the quench amplitude is very small compared to its width, the Ward identity (4.6) 
implies that the system’s energy will be very nearly conserved, and thus we anticipate 
final configurations that are very close to the initial state. In particular, an initial state 
characterized by a small black hole can remain within the small black hole branch. 

An example of this kind of quench is shown in figure 12. There the perturbation 
width is taken to be the same as in the example of previous section, f = 0.168, but 
now the amplitude is chosen to be 5 = 0.02. 
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Figure 13. A slow quench with large amplitude: S — 0.5 and f = 1.19. 

With these parameters, the area of the event horizon in the final state is only about 
2.5% larger than that of the initial state. Pictorially, on the scale of figure 7 this quench 
results in the red dot which marks the location of the initial state moving imperceptibly 
to the left after the perturbation. At late times, the ring-down of (O) and (T^x) is far 
more pronounced in this case, with easily identifiable oscillations shown in figure 12. 
The late time behavior of this quench is the focus of figure 8, which highlights the 
role of the dominant quasi-normal mode. The late time behavior of (O) displays very 
similar behavior. In line with expectations, unlike the large amplitude quenches this 
class of perturbation is sensitive to the particulars of the initial state. For example, 
for the quench shown in figure 12, we find that Ttherm ~ 2 in units of /o—but the 
same perturbation in an initial state with about 13.5 times the energy density has a 
thermalization time Ttherm 0.3, more than 6.6 times faster. 

Of course there are a variety of perturbations one can perform which result in a 
final state thermodynamically similar to the initial state. As mentioned above, these 
all share the feature that the perturbation amplitude is relatively small compared to 
its width. We will see in detail below how to quantify relatively small. At present, we 
turn our attention to a similar example of thermalization along a branch. This is the 
large amplitude, slow quench shown in figure 13. In this case, (5 = 0.5 and f = 1.19, so 
that the ratio of the amplitude to width here is about 3.5 times the previous example. 
Accordingly, we anticipate a greater change in energy density and that this perturbation 
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Figure 14. The final state energy density for various perturbations with either 6fo or r 
fixed. All quantities in this figure (and those that follow) are measured in units of /o 


will thermalize more quickly. From the late time behavior we find that this is indeed 
the case, with Ttherm ^ 0.86. 

Our code is able to evolve even slower perturbations, including those whose widths 
are many times larger than their amplitude. Predictably, in this case the various one 
point functions respond by being gradually deformed from (and then nearly returned 
to) their initial vales. For these adiabatic perturbations, the “bumpy” features in the 
late time behavior shown in figure 13 disappear, and the exponential decay at late times 
is difficult to detect due to numerical accuracy. This can be viewed as a consequence 
of the fact that for such slow perturbations, the quench width is much larger than the 
time scale defined by the final state’s lowest lying quasi-normal mode. 


7 Perturbation Analysis and Dynamical Phase Structure 

7.1 Parameter Dependence and Scaling Regimes 

To better understand the featnres of different qnenches, we look in more detail at the 
dependence of the final state on the perturbation parameters. In figure 14, we show 
several results for the dependence of the final state energy density on 5 for fixed f and 
vice-versa. In figure 14(a), we compare perturbations with many different amplitudes 
but one of three fixed qnench widths: f = 0.168,0.433 and 0.838. In fignre 14(b), we 
fix the pertnrbation amplitnde at h = 0.05, 0.5 or 1 and consider quenches with many 
different widths. 

The fixed r scenario is studied closely in figure 15. For sufficiently small S and 
f, the change in (Tu) across the quench is well approximated by a quadratic in the 
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5fo .^fo 

(a) T = 0.168 fixed (b) r = 0.838 fixed 

Figure 15. Analysis of two fixed r cases. In the small 5fo limit, both cases can be well ht 
with quadratic functions of 5fo (insets). In 15(a), the large SJq region favors a power law ht 
with (Tit) ~ Around 6fo — 0.7 in 15(b), a linear behavior (Tit) ~ <^/o appears, while 

some deviations are seen around (5/o ^ 1 in this slow quench. 

amplitude of the perturbation. In fact, this behavior manifests over a wide range of 
quench amplitudes providing that the quench is very fast, so that f is very small (or 
r <C /o)- As the quench width is increased, however, deviations from this simple scaling 
become more pronounced. In hgure 15(b), in which f is no longer small, the quadratic 
amplitude dependence of the final state energy density is only present for 5 < 0.1. 
Outside of this small amplitude region we find that in this case (T^) ^ Sfo. 

In figure 16 we analyze the dependence of the final state energy density on the 
quench speed with fixed S. Again, there is a pronounced scaling regime evident when 
the quench is fast. In this case we find that (Tu) asymptotically behaves as l/r^ even 
when d/o is comparable in size to /q. The scaling behavior degrades as f becomes 
large, as is obvious from the lower left corners of the plots in figure 16. This suggests 
that it is roughly the duration of the quench as compared to the characteristic scale 
of the theory /o that determines whether we are in this scaling regime, as opposed to 
the size of the width relative to the quench amplitude. In figure 16(b) we consider 
perturbations with fixed small amplitude. As shown in the inset, in this case there 
appears to be an intermediate region around f ~ 1 where the {Tu) ^ scaling 

momentarily returns. However, as the quench width is further reduced we find that 
this scaling is not maintained, and the functional form appears to approach {Tu) r\j 1/t 
in the adiabatic limit. This adiabatic behavior at small 6 is consistent with what was 
found in the absence of a confining potential in [14]. 

Combing the lessons of this subsection’s results, we find that {Tu) oc ((5 /o/t)^ when 
the quench width is small compared to all other scales. This implies that a constant 
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(a) 6fo = 1 fixed 


(b) 6fo = 0.05 fixed 


Figure 16. Analysis of two fixed r cases. In both cases, the small r asymptotic behavior 
is given by (Tu) 1/r^. In 16(b), there is a distinct intermediate region near r 1 where 
(Ttt) ~ 1/r^ once again before transitioning to {Tu) ~ 1/r in the adiabatic limit (inset). 


(fu) line on the quench parameter plane is given by Sfo ~ r. The scaling appears to 
be approximately respected even for small amplitude perturbations, which is consistent 
with the linear behavior found in the bottom left corner of the dynamical phase diagram 
in figure 17. Perhaps not surprisingly, this behavior is also consistent with a result first 
observed numerically under fairly different circumstances [14]. It was subsequently 
demonstrated [16-18] that such a scaling is universal in the sense that it only depends 
on the near boundary features of the bulk theory, which is asymptotically AdS in many 
applications relevant for holography. Roughly, the idea is that very abrupt processes 
in the UV do not probe very deeply into the bulk, as the quench concludes before the 
disturbance has time to propagate far from the boundary. Since all the interesting 
features of our model appear deep in the IR, they play no role in the dynamics of 
this fast quench regime and it is therefore natural for our model to exhibit the same 
universal scaling behavior. 

7.2 Dynamical Phase Structure 

The coarse features of the quench dynamics can be efficiently summarized in a dynam¬ 
ical phase diagram, which we present in figure 17. To construct this diagram, we begin 
in the smallest black hole with which we can reliably evolve arbitrary perturbations, 
and scan the (f, 5) parameter plane. For each of the many perturbations, we follow the 
evolution until a static final state solution is obtained, and then use figure 5 to decide 
if that solution describes a large or small black hole. 

The resulting phase diagram has several noteworthy features. First, the thick 
dashed line separating the perturbations which result in a small black hole from those 
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Figure 17. Dynamical phase diagram for the small black brane initial state Xh/^c = 3.23. 
The dashed line (lower) marks the location of a crossover between perturbations which result 
in a small black brane and those which produce a large black brane final state. The dotted 
line (upper) is the location above which the large black hole is thermodynamically favored 
over the thermal gas solution. Both lines behave like 5 ^ f for f < 0.1. 


that end up in a large black hole likely marks the location of a crossover, as we have not 
been able to find any non-analytic behavior across this transition. To look for the sort 
of non-analytic behavior we have in mind, one can consider the set of susceptibilities 
that describe the system’s dynamical response to various perturbations. For example, 
derivatives of the form 






(7.1) 


where AS is the change in the system’s energy density and the subscript is an instruc¬ 
tion to take derivatives along the direction where that parameter is held hxed. 

In our quenches, these susceptibilities appear to be continuous and smooth as one 
approaches the transition line. This can be inferred from the plots shown in figure 18 
which shows an example typical of what we find when we differentiate the final state 
energy density with respect to the quench parameters at fixed 5 ov f. In general we 
expect that most dynamical observables in our model such as thermalization times and 
entropy production depend smoothly on AS = (T**) final — (T’^Oinit, and thus we do 
not anticipate critical behavior to arise in other observables in the vicinity of this line. 
This is analogous to what happens in weak field perturbations of global AdS, in which 
the line separating small and large black hole final states was also consistent with a 
crossover [8]. 

Finally, the functional form of this dashed curve appears to be suitably described by 
a simple power law in the “small-fast” regime. Here the data can be well fit to a function 
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Figure 18. Derivatives of {Tu) in the final state with respect to the quench parameters 
are smooth across the line separating large and small black brane final states. The vertical 
dashed lines mark the location of the “critical” quench parameter which brings the initial bulk 
solution to the minimum temperature black brane. Derivatives taken at other fixed values of 
6 and f show similar behavior. All quantities are measured in units of /q. 

of the form 5 ^ f, which is in line with general expectations as we discussed above. 
Independent of its precise functional form, we find that the crossover line monotonically 
increases in 5 as f increases for the parameter space covered by our study. 

8 Discussion 

The line of research initiated in the present work provides another step forward in the 
holographic study of gauge theories that are qualitatively similar to SU (3) Yang-Mills 
in 3-1-1 dimensions. Our gravitational model is characterized by a non-trivial (and care¬ 
fully tuned) dilaton potential which introduces an additional level of complexity to the 
gravitational infall calculations dual to boundary theory thermalization. Specifically, 
the gravitational solutions of our model all involve a running dilaton in the bulk whose 
profile encodes the breaking of conformal invariance in the dual gauge theory. 

The fact that these “initial state” static solutions have nontrivial dilaton profiles 
requires a careful treatment of their near-boundary behavior in order to achieve stable 
evolution. We have applied a numerical technique based on finite difference integration 
to successfully evolve a broad class of perturbations, solving the fully backreacted 
Einstein-dilaton equations for the duration of the perturbation/equilibration process. 

By studying the dependence of the final state on the parameters of the quench, we 
have collected several interesting lessons which we expect to prove insensitive to the 
precise details of the dilaton potential. Among these are the observation of a rapid 
transition from the perturbation (initial condition) dominated regime to the quasi- 
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normal mode dominated regime of the quench dynamics. This behavior appears to be 
ubiquitous in holographic realizations of thermal quenches, and likely reflects the fact 
that for solutions with sufficiently large black hole horizons the thermalization time is 
governed by the Hawking temperature Ttherm ~ l/T (which is now the dominant length 
scale). The extent to which this continues to be true in the extremal limit of our model 
is largely an open question. From figure 9 it is clear that the decay width of the lowest 
lying quasi-normal mode has already departed from the simple conformal temperature 
scaling for final states with temperature near T^. However, as we illustrate in section 
6 .2, our perturbations which remain on the small black hole branch also manifest this 
quick transition to the linear regime. An interesting question (intimately related to the 
second question concluding section 6) is then where might one expect to see deviations 
from this “ubiquitous” behavior? 

As we move further along the small black hole branch, the length scale introduced 
by the presence of the dimensionful source becomes increasingly important. This sug¬ 
gests that a sensible guess for where deviations from the rapid transition to the linear 
regime appear is for those perturbations whose final state has an energy density which 
is small compared to the source—in other words {Tu )/^ 1- The smallest black 
hole which we perturb in the present work has {Tu )/~ and thus it is perhaps 

not surprising that we see no evidence for a new non-linear regime sensitive to the 
presence of the confining potential. Pushing the reach of our numerics closer to the 
extremal solution is an ongoing direction of our research which we hope to report on 
in the future. 

Ultimately, we would like to go beyond the small black hole limit we have studied 
in this work and consider perturbations of the horizon-less zero temperature solution 
directly. The most compelling motivation for this is to determine whether or not the 
diverging dilaton potential characteristic of a broad class of holographic models of QCD 
is sufficiently repulsive to give rise to scattering type solutions which never result in 
black brane formation. If this is indeed the case, our dynamical phase diagram would 
gain a new line dividing those perturbations which result in black brane formation from 
those which do not. Finding the associated scattering solutions would have interesting 
implications for the dual gauge theory, suggesting a class of perturbations that the 
strongly coupled matter can not thermalize. 

Moreover, the boundary of these “unthermalized” perturbations in the dynamical 
phase diagram would be interesting in its own right. By analogy with more familiar 
examples in asymptotically flat space, one might hope to find critical behavior akin to 
Choptuik phenomena [39], in which a bulk scaling solution appears on the boundary 
of the perturbations that do and do not form a small black hole. Solutions of this 
sort are typically accompanied by various power law scalings characteristic of a second 
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order phase transition. In our model, for example, this could manifest as the final state 
energy density assuming the form Af |f ~ (5 — where Sc is the amplitude of the 
critical perturbation and 7 is a critical exponent quantifying the universality class of 
our model. The exponent may be different from that of Choptuik, as the small black 
holes in our theories (unlike AdS) depart importantly from flat space black holes. 

The appearance of a universal scaling regime in the fast quench limit is another 
noteworthy output from our calculations. As previously discussed, this scaling regime 
was anticipated on very general grounds, and its manifestation in our abrupt quench 
data is in some ways an encouraging check on our numerics. It is interesting to wonder 
how this scaling might be effected by the potential barrier inherent to the zero tem¬ 
perature solution. In so far as the universality of this scaling depends only on the UV 
features of the bulk solution, it is likely that quenches of the extremal solution whose 
width is much smaller than the scalar source will again result in a final state whose 
energy density satisfies (1.1). However, because the arguments for this scaling behavior 
are closely tied to the early time dynamics of the quench, it remains unclear whether 
the increase in energy density will manifest in black brane formation or a non-thermal 
scattering solution. 

Investigating the properties of other probes which can be used to characterize our 
perturbations is another interesting future direction. In [15, 40] a variety of non-local 
probes were used to measure the approach to thermal equilibrium. These include 
various two point correlation functions, Wilson loops, and the system’s entanglement 
entropy. The advantage of these non-local probes is that they are capable of moving 
beyond the binary thermal/non-thermal characterization of the perturbation, as they 
are sensitive to the scale dependence of the thermalization process. In the gravitational 
picture, they accomplish this by sampling the geometry away from the UV boundary. 

To understand why this is true, it is instructive to consider the equal time two-point 
correlator for some gauge invariant boundary operator with large conformal dimension. 
In the semi-classical limit, this correlation function can be computed holographically 
by calculating the length of the bulk geodesic that connects two spatially separated 
points on the boundary. As the separation distance between these points increases, the 
bulk geodesic droops increasingly deeper into the IR. By measuring the deviation of 
the length of this bulk geodesic from the value one obtains in the final thermal static 
state as a function of boundary separation and time, one arrives at a picture of the 
approach to thermal equilibrium at different length scales. Constructing this picture 
in our model is currently in progress. 

Following the discussion above we may put forth the following qualitative expec¬ 
tations, as a function of the bulk holographic theory: 
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• For theories that are simple (and smooth) RG flows between two CFTs (or hy- 
perscaling violating geometries), the physics is characterized by a single mass 
scale. In such cases, the T = 0 theory is gapless and at T > 0, the theory is in 
the black hole phase. There is no phase transition or very fast crossover. The 
typical diagram that maps the black-holes of this theory is as shown in the left 
flgure 19 where the temperature of the black hole solution is plotted against the 
value of the driving bulk scalar at the horizon, (ph. In such cases, generically the 
thermalization time is expected to be proportional to the temperature T. 

• For theories that are RG flows to a gapped IR theory as the case studied in this 
paper, the situation is different. Such theories are confining^ and have a first- 
order phase transition at T = Tc ~ A to the deconfined plasma phase (A is the 
confinement scale of the T = 0 theory). The typical diagram that maps the black- 
holes of this theory is as shown in the central flgure 19 where the temperature 
of the black hole solution is plotted against the value of the driving scalar at the 
horizon, ph. There is a large stable black-hole branch to the left and a small 
unstable black-hole branch to the right. 

In the confined phase, T <Tc, the spectrum is discrete and real to leading order 
in l/Nc- The imaginary parts of various quasinormal modes are of order 1/Nc, 
and therefore the thermalization time, is of order A“^A^c if the energy density 
injected into the system is much smaller than T^. In particular, using the tree- 
level bulk equations of motion we do not expect the system to thermalize as the 
imaginary parts are zero to that order. We may therefore expect a Ghoptuik-like 
phase transition in that case. 

In the black hole phase and for energies near the phase transition the thermaliza¬ 
tion time is expected to be of order Tc ~ A. Finally for energy densities ^ A^ 
we expect the thermalization time to be set by the (final) temperature T. 

• There are intermediate cases in which the T = 0 theory is gapless in the IR, but 
with a fast crossover or phase transition to the UV regime. Such theories are 
always deconfined (according to the Wilson loop test) at T = 0 [41, 42]. 

The typical diagram that maps the black-holes of this theory is as shown in the 
right flgure 19 where the temperature of the black hole solution is again plotted 
against the value of the driving scalar at the horizon, ph- There is a large stable 
black-hole branch to the left and a small stable black-hole branch to the right, 
while there is also an unstable black-hole branch in the middle. In such a theory 

^In Einstein dilaton theory a confining theory is always gapped and vice-versa [41, 42]. 
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Figure 19. Left: Typical {(f>h,T) diagram of the gapless/non-confining theory. Center: 
Typical {4>h,T) diagram of a confining and gapped theory. Right: Typical {4>h,T) diagram 
of a non-confining and gapless theory with a first order phase transition. 

there is a (continuous) phase transition at T = O'*' to the small black hole phase, 
followed by a first order phase transition (or a fast crossover) to the large black 
hole phase®. 

For very small or very large final state temperatures, it is the temperature that 
sets the thermalization time, but in the intermediate region the characteristic 
scales of the theory enter into the thermalization time. Moreover, for small energy 
density we do not expect a Choptuik-like threshold in this case. 

An important question is the implications of our results for heavy-ion collisions. 
The fast thermalization of sufficient high-density initial states should be comparable 
to the processes studied here upon translating appropriately the quench parameters. 
Of course it should be kept in mind that the heavy ion collisions are anisotropic in 
space, but there are good reason to believe that this is not very important for initial 
thermalization but more important for the subsequent evolution of the plasma. In this 
respect our numerical results on thermalization time should provide reliable estimates 
for the analogous heavy-ion thermalization process. Most importantly, the part of the 
physics that is not described here, namely the boundary between thermalization and 
non-thermalization will also provide important clues for thermalization in the less un¬ 
derstood pp collisions, where recently CMS reported the first ever evidence for collective 
effects, [44, 45]. 

Note Added 

Immediately before this work was submitted, several papers [46-48] focusing on ther- 
malization processes in nonconformal theories/backgrounds simultaneously appeared 

®Such RG flows have two scales at T = 0 and therefore a dimensionless parameter. An example is 
YM with quarks of mass m, and such a black hole diagram was found for example in [43] (see figure 
22 of that paper). 
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on the arXiv. These works have some overlap with our results on thermalization times 
in a non-conformal field theory. 

Specifically, [46] focuses on the important role played by the low-lying quasi-normal 
modes in the approach to thermal equilibrium. The fluctuation problem they consider 
is a linearized fluctuation problem in backgrounds like Af = 2* which are gapless, with 
hyperscaling violating asymptotics in the IR. The authors emphasize the fact that 
even in states of non-conformal field theories with large deviations from conformality, 
the lowest lying quasi-normal mode approximately scales linearly with temperature. 
This is certainly true in many examples of holographic matter in the high temperature 
plasma phase, and in the corresponding states of our model as well (figure 9 ). However 
this approximate scaling is almost certainly violated in more phenomenologically viable 
models of the QGP, where the interaction measure is strongly peaked near Tg. Indeed 
this scaling is also violated noticeably in the phenomenological models investigated in 
[48]. 

We believe that the subsequent claim made by the authors of [46], that “the ther¬ 
malization time is generically set by the temperature, irrespective of any other scales, 
in strongly coupled gauge theories” is too strong and valid only in theories with a single 
dynamical scale which are non-confining and gapless at T = 0. Such systems have a 
phase transition to a black hole phase as soon as T > 0. For reasons discussed in [11] 
and in the body of our text, we find it plausible that the scale introduced by the mass 
gap in our model introduce a new dynamical regime which is distinct from the QNM 
ringdown. Therefore, even the assumption made by those authors that thermalization 
times are approximately bounded from above by the lowest lying quasi-normal mode is 
called into question. Evaluating these claims necessarily requires going beyond the lin¬ 
earized gravitational equations, which is the approach we have adopted in the present 
work. 

On the other hand, the authors of [48] focus again on linearized quasinormal modes 
of a minimally coupled scalar. The difference now is that the background is an Einstein- 
dilaton theory with a fast cross-over (and in one example a phase transition). The main 
differences from our theory is that their theories are gapless at zero temperature while 
ours are gapped, and in our case the relevant scalar is the same that participates in the 
vacuum solution. Some of their phenomenological formulae (like the connection of the 
thermalization time to the speed of sound) do not work well below Tc in our case, and 
the reason may be the differences stated above. 

Einally, the work of [47] is closest in spirit to the present work, as these authors 
perturb their system and they follow the non-linear evolution. The main difference lies 
on the system to be perturbed. In their case the system is A/" = 4 plasma at finite 
charge density or in the presence of an external magnetic field. 
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APPENDIX 


A Holographic Renormalization and Boundary Theory Cor¬ 
relators 

The near boundary analysis of section 2 . 2.1 is convenient for developing an algorithm 
for solving the Einstein equations, but the in-going null coordinates are ill suited to con¬ 
structing the generating functional of the dual gauge theory. Following [32] , the analysis 
of the on shell action near the boundary is most directly performed in Fefferman- 
Graham coordinates, in which the radial direction is orthogonal to the boundary direc¬ 
tions. 

In Fefferman-Graham coordinates, the metric takes the form 

ds^ = ^ ^ Qij dxMa:-^, (A.l) 

and the metric and scalar permit the following expansions: 

9{t,p) [g{n){t) + h(n){t)\og p + \n){t)\og^ p + ■ ■ ■] p'^, 

n=0 

^{t,p)='^[‘P(n){t)+'lp{n){t)logp + 'lpn{t)log^P+---] (A.2) 

n=0 

Substituting these expansions into the Einstein equations and solving them order by or¬ 
der in p allows one to determine many of the coefficients in the expansions algebraically. 
The primary exceptions are the leading coefficients of the normalizable modes, 5 ^( 4 ) and 
(p){ 2 ) which can only be determined given the full radial profile in the bulk. Neverthe¬ 
less, the near boundary analysis does constrain these undetermined coefficients, a fact 
realized in the boundary gauge theory by the existence of Ward identities. 

The regularized on shell action ( 2 . 1 ) can be written 

<5r = -^/ dpdVC^T( 9 ?) - ^^dVC= 7 /C, (A.3) 

where Einstein’s equations have been used to eliminate the Ricci scalar. This action 
exhibits the following divergences in the limit e ^ 0 : 

~ g'^(0)'^b(0)] + 2A(o)noV^(o)) + 0(e°) 


(A.4) 




R 


2«:2 


d^a;y-^(o) 
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To construct the appropriate counter terms, the regulated action must be expressed 
in terms of the helds living on the surface p = e. These helds are the pullback of the 
metric, 7 (t,e), and the scalar (p{t,e). Performing this inversion yields 




fp=e 


6 + -R[y] + + log e ( F 4 + -(p^R[^] 


9 


-\{^jRi,b]R'=b] - +^[7, v] 


(A^5) 


where A contains 0(e°) finite counter terms, and the coefficient F 4 depends on the 
details of the higher order terms in the scalar potential. For the potential in (2.6), it 
reads 




(A.6) 


In the absence of an organizing principle such as supersymmetry, the finite counter 
terms are left unfixed and thus lead to scheme dependent ambiguities in the correlation 
functions. 

Once the counterterms have been identified, one can form the subtracted action 

like 

'S’sub = <5r + Sc, (A.7) 


and the renormalized correlation functions are then computed as follows: 


(O) =hm 

' ' e-i^O 

(T,,) =hin 


(I 1 ^>gsub \ 

W \/^ / ’ 



(A.8) 

(A.9) 


where Tjjfy] is the stress tensor of the theory at p = e. This boundary stress tensor is 
generically the sum of the contribution from the regularized action, and the contribution 
due to the presence of the counterterms. The regularized action gives 

Fb\ = -h - I<1») = I b.7« - (A.IO) 

Z 

with Kij the extrinsic curvature of the regulating surface and K = its trace. To 

obtain the contribution from the connterterms, it is convenient to first catalogue the 
metric variations of a boundary action of the form 

s^ = J d^x^(^A + BR[^]+CR^[^] + DRabb]R'^'’b] + EpDY, (ATI) 
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where A, B, C, D and E are arbitrary scalar functionals independent of the metric. 
Obviously this action contains (A.5) as a special case. With a bit of effort, one can 
show that the metric variations of yield 

2 

Ttf [ 7 ] = =lij {A + BR + CR^ + DRa,R^^ - V,(^<^)VV) 

- 2BRij + 2 VjVji? — 2'^ij 

- ACRRij + AViVj{CR) - A-iij V^{CR) 

- ADRi ^Rja + AVkV^i{DRj) - 2V\DRij) - 2^,^ 

+ 2V(i(£'(/p)Vj)9?. (A. 12 ) 

All contractions and curvatures in 7)^ implicitly refer to the metric 7 . From the terms 
in (A. 12 ) with coefficients determined by comparison to (A.5), it is straightforward to 
obtain the contribution to the boundary stress tensor resulting from the counterterms. 

Performing this maneuver, summing the result with (A. 10) and then inserting 
the on-shell near boundary expansions from (A. 2 ) into (A. 8 , A. 9), one obtains the 
renormalized one point functions: 

2 / \ /16 

= 2 (^4 (yC( 2 ) + (/3(o) j ^ J + (A.13) 

1 / \ 1 1/16 

K^{Ttt) =2g^4)tt - g (^8(/P(2) - (/(o)jv?(o) - g</(o)^ + ^ J + 

(A.14) 

2 1 / \ 1 11 /208 \ 

K {T^t) = g9(4)« + Jgp<9(2) + 6.^(0)j(/3(0) + g<9(0) + I ^ j + 9A. 

(A,15) 

The schematic notation dA refers to the contributions coming from the set of finite 
counter terms contained in A. These contributions are scheme dependent, and will 
henceforth be neglected for simplicity. Note that in deriving these expressions the 
boundary metric has been assumed to be flat, g{Q) = gij. Written in terms of the 
near boundary expansion coefficients, it is straightforward to demonstrate that these 
one-point functions respect the anticipated Ward identities. For example (in units of 

2 ( \ 2 1/16 
= 3 (4<yi5(2) + </( 0 )j‘y^( 0 ) + ^ ^ 

2 1 /16 

= V^(o) {O) + ^j (A.16) 
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demonstrates the breaking of conformal symmetry in the presence of a dimensionful 
source in terms of the classical result (hrst term) and terms due to the matter anomaly 
in four dimensions. Similarly 

2 / \ /16 

V*(rtt) = -(^4(^(2) + (p(o)j<P(o) + -V^(o) = ^(0){O) (A.17) 

describes the change in the system’s energy in terms of the work done on it by a time 
dependent source. The derivation of these identities requires additional constraints 
which are easily obtained from the equations of motion. They relate, for example, time 
derivatives of the undetermined near boundary coefficients. 

Because the numerical computations directly access the boundary expansion coef¬ 
ficients given in (2.20), it is convenient to relate these coefficients to those appearing 
in (A.2). The coordinate change is given by 

dx^ dx^' 

which implies two particularly useful equations: 

\ ^ - v'‘^A - ^v'z\ (A.19) 

0 — vv' A —^ {yz! + iu'). (A.20) 

These equations can be solved perturbatively in p to obtain z{t,p) and v{t,p). One 
straightforward method to this end is to expand the ingoing null coordinates in powers 
of p, like 


Z{t,p) =p + '^p"'{Sn{t) + Sn{t) log p), 

n=2 

v(t,p) =t + ^p^{cn{t) +Cn{t) logp), 


(A.21) 

(A.22) 


n=l 


and substitute these expansions into (A. 19) and (A.20). The resulting system can be 
solved order by order in p for the coefficients Snit) and c„(t). This procedure yields 


p)=p- ^/o + ^/o/o ^ 

- ^«4p^logp + 0(p®), 


32 


16 

404 + «4 + ^/o-^ 


u(t, p)=t- p - ^fl p^ + ^/o/o + 0{p^), 


(A.23) 

(A.24) 
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and upon inserting these expansions into (2.18) and ( 2 . 20 ) and regrouping terms, one 
directly obtains the independent coefficients ip(Q),ip( 2 ) and g{ 4 )u in terms of /o,/ 2 , and 
04 . The result is 



(A.25) 


which can be used to write the one-point functions in terms of the coefficients obtained 
directly from the numerical routines, as in (4. 2-4. 4). 

B Initial Data and Convergence 

The initial data we wish to perturb and evolve in time are the solutions to the static 


equations of motion given by setting n-derivatives to zero in (2.8-2.12). These solutions 


are constructed by integrating the static equations from the horizon to the boundary 
and matching to the boundary behavior of solutions in the gauge (^{y) — 0. This com¬ 
putation can be performed at very high precisions without much effort. The results are 
then exported on the discretized grids desired for time evolution. Since the numerical 
codes used to produce and evolve the static solutions are distinct, it is an important 
check on our numerical package that the unperturbed initial state can be stably evolved. 
This procedure also helps us determine suitable grid sizes for satisfactorily suppressing 
numerical errors. 

When the energy density of the initial state background is large, the dilaton is 
small throughout the bulk and it is fairly easy to achieve robust evolution. These 
static solutions can be evolved with only moderately small grid sizes, and numerical 
errors are very small. The scaling of the discretization error is second order in /S.z, as 
desired. 

As the background scalar becomes larger, however, we find that even simply main¬ 
taining the static initial state becomes a challenge. In particular, the difficulty dras¬ 
tically increases as the black hole size decreases as it heads towards the small black 
hole branch. In figure 20, time evolution of a static solution with Xr/K = 3.23 is 
shown. This corresponds to almost the smallest black brane in which we can reliably 
perform many distinct perturbations. We divide the domain between the boundary 
and the initial black hole horizon Zif|stat into N intervals, and hence the grid size is 
Az = Znlstat/A. In the left panel is the time evolution of {O). When N = 200 and 
400, there are some numerical oscillations which persist on the order of Az, which is 
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(a) Static (O). (b) Log plot of static S{Ttt)- 

Figure 20. Convergence of static time evolution by changing the grid size for an initial 
data set with Xh/^c — 3.23. In 20(b), the magnitude of the difference of {Ttt{v)) from the 
input value is shown, S{Tu) = \{Ttt{v)) — (T)t(0))|/|(T)t(0))|. The blue, purple and green 
lines correspond to N = 200, 400, 800, respectively. In the late time, the green lines do not 
oscillate and converge to the static value exponentially. 

larger than the desired numerical error. When the grid size is decreased by one half, 
although oscillations remain, the decrease of their magnitude is demonstrably second 
order in Az. This implies that the presence of these oscillations is not indicating any 
real numerical instability, but is instead noise. We observe that this noise is largely 
insensitive to the size of Av. For the initial data shown in figure 20, we find that 
the noise is heavily suppressed and the late time behaviors of both panels converge 
to constant values when N = 800. These results suggest that for large background 
scalars, taking a very small grid size is necessary in order to rid the computation of 
unwanted noise. We have also made several tests with larger Xh/Xc, and in those cases 
the required smallness of the grid size for suppressing the noise continues to quickly 
increase. The static evolution of (T^x) behaves similarly to (O). 

These oscillations are far less pronounced in {Tu), as indicated by (4.6). In the right 
panel of figure 20, we plot the magnitude of the difference of (Tu) for the unperturbed 
solution evolved in time from the input value. Notice the difference of scale compared 
to the left panel. In this logarithmic plot, quadratic convergence against the grid size 
of the finite discretization is evident. We also checked the bulk constraint equation 
(2.12), and verified that this too converges quadratically. 

One consequence of evolving our unperturbed initial state solutions in time is that 
any numerical irregularities in the solution are damped by the presence of the horizon 
as the numerical solution “rings down” to a numerically stable configuration. This 
procedure thus cleans the numerical data for static initial states from the discrepancy 
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between using different numerical codes for static and dynamical computations. In 
practice, before we perform any sort of quench, we typically allow the numerical data 
for the initial state to evolve unperturbed for a bit to clean the noise, and then compute 
the evolution of the perturbed state with suitably small (quench dependent) grid sizes. 

C Geometrical Aside 

C.l The Apparent Horizon 

The formulation of our gravitational problem relies crucially on the notion of an ap¬ 
parent horizon. This horizon is primarily important because unlike its more familiar 
sibling the event horizon, it is not teleological. This is to say that the location of the 
apparent horizon can be determined on each time slice, whereas the event horizon can 
only be deduced once the final state of the geometry is known. 

An operational definition of an apparent horizon is a spacelike surface on which an 
outgoing null congruence normal to the surface has zero expansion. Following [49] one 
may study this expansion 9 via a null vector tangent to an outgoing null geodesic, k. 
In the basis provided by the coordinates of (2.7), such a vector is given by 

k = k^d,^d,-^Ad,. (C.l) 

This vector is not affinely parametrized, which means that it satisfies the geodesic 
equation 

k^V ^k,y — nki, (C.2) 

for non-zero k. It is easy to show that in the present case, 

(0.3) 

where the prime denotes differentiation with respect to . 2 . Because k is non-zero, the 
expansion equation is modified as follows: 

e = e^ (V^F - k) . (C.4) 

In this expression, the exponential pre-factor is necessary to convert to an affine 
parametrization. As it is manifestly positive and non-zero, it will play no role in 
the present discussion. 

From (2.7), (C.l) and (C.3) a short calculation reveals that in this background 
ansatz ^ 

e = dj:-^Adj:, (c.5) 
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and combining this result with the definition of the outwards directed derivative d_|_ 
from (2.13), and the requirement that the expansion vanish at the location of the 
apparent horizon zh leads to this section’s main result: 

d+SL„ = 0. (C.6) 

Apparent horizons have several interesting properties that make them particularly well 
suited to the studying of dynamical processes in gravity. Specifically, for static space- 
times with black holes in the bulk, the apparent, event, and Killing horizons all coin¬ 
cide. Moreover, it is possible to show on general grounds that in the non-static case 
an apparent horizon will always he within the event horizon. This suggests that an 
apparent horizon provides an IR cutoff for the numeric problem that is both natural 
and practical, as well as available on each time slice. 


C.2 The Event Horizon 


When the dynamics is such that the gravitational system settles back into a static 
state, as is the case in the holographic thermalization processes we study here, it is 
also straightforward to compute the location and area of the event horizon over time. 
For this one needs only the observation (mentioned above) that the event and apparent 
horizons coincide in static spacetimes, and that the event horizon necessarily travels 
along a geodesic in the bulk. 

The line element traversed by a light ray in the ingoing null coordinates satisfies 

Q = -Adv^-^Mz, (C.7) 

z 

which implies the following geodesic equation for the location of the event horizon, 


^EH — — 2^EH^(^EHW)) 


(C.8) 


to be solved subject to the boundary condition z-^niv oo) = zxr{v oo), where 2 :ah 
is the location of the apparent horizon. In practice, equation (C.8) can be integrated 
backwards in time along the geodesic using the numerically determined metric function 
A at each step. 
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